/* ====================================================
#   Copyright (C) 2021-2025  All rights reserved.
#
#   Author        : Jinjing Li
#   Email         : josiahleeoaa@outlook.com
#   File Name     : libAlphaN.C
#   Last Modified : 2025-05-09 00:13
#   Describe      : Implementation of (alpha,n) reaction neutron generator
#                   - Updated parameter lists (removed unused params)
#                   - Added detailed function documentation
#                   - Improved code comments
#
#   Features:
#   - Processes NEUCBOT nuclear reaction code outputs
#   - Generates neutron and gamma spectra
#   - Supports multiple excited states
#   - Provides angular distribution sampling
#
# ====================================================*/

#include "libAlphaN.h"
#include <numeric>
#include <sys/stat.h>
#include <cerrno> // for errno

using namespace std;
using namespace CLHEP;

// ==============================================
// Main Generator Functions (Top Level Interface)
// ==============================================

/**
 * @brief Generate (alpha,n) reaction events for single excited state
 */
void generator_alphan_basedon_sundean(int isotope, const string &target, int state_idx,
                                      unsigned int nEvents, FILE *stream,
                                      const string &database)
{
    string keyname = get_alphan_key(isotope, target, state_idx, database);
    vector<double> neutron_energy_vals;
    // Sampling neutron energy from root files
    get_neutron_energy_from_sundean(keyname, nEvents, neutron_energy_vals);

    TRandom3 genor(0);
    gRandom->SetSeed(0);

    fprintf(stream, "# File generated by u_generators.\n");

    vector<string> key_mode;
    vector<double> intensity, energy;

    if (state_idx != -1)
        get_gamma_energy(target, state_idx, key_mode, intensity, energy);

    double cosTheta1, theta1, azimuth1;
    for (unsigned int i = 0; i < nEvents; i++)
    {
        cosTheta1 = genor.Uniform(-1., 1.);
        azimuth1 = genor.Uniform(2.0 * M_PI);
        theta1 = acos(cosTheta1);

        double Ek = neutron_energy_vals[i];
        double momentum1 = sqrt(pow(Ek / GeV + neutron_mass_c2 / GeV, 2) - pow(neutron_mass_c2 / GeV, 2));
        dirs p1(momentum1, theta1, azimuth1);

        if (state_idx == -1)
        {
            fprintf(stream, "1\n");
            fprintf(stream, "1\t2112 0 0 %e %e %e %e 0\n",
                    p1.x_norm, p1.y_norm, p1.z_norm, neutron_mass_c2 / GeV);
        }
        else
        {
            vector<double> Eout;
            generate_deexcitation_gamma(key_mode, intensity, energy, Eout);
            int num_g = 1 + Eout.size();
            fprintf(stream, "%d\n", num_g);
            fprintf(stream, "1\t2112 0 0 %e %e %e %e 0\n",
                    p1.x_norm, p1.y_norm, p1.z_norm, neutron_mass_c2 / GeV);

            for (int ig = 0; ig < num_g - 1; ig++)
            {
                cosTheta1 = genor.Uniform(-1., 1.);
                azimuth1 = genor.Uniform(2.0 * M_PI);
                theta1 = acos(cosTheta1);

                double Eg = Eout[ig];
                double momentum2 = Eg / GeV;
                dirs p2(momentum2, theta1, azimuth1);
                fprintf(stream, "1\t22 0 0 %e %e %e %e 0\n",
                        p2.x_norm, p2.y_norm, p2.z_norm, 0.);
            }
        }
    }
}

/**
 * @brief Generate (alpha,n) reaction events for multiple excited states
 */
void generator_alphan_basedon_sundean(int isotope, const string &target,
                                      unsigned int nEvents, FILE *stream,
                                      const string &database)
{

    vector<double> neutron_energy_vals[MAX_LEVELS];
    vector<string> key_mode[MAX_LEVELS];
    vector<double> intensity[MAX_LEVELS];
    vector<double> energy[MAX_LEVELS];
    vector<double> branch_ratios;

    get_gamma_energy(target, isotope, database, key_mode, intensity, energy, branch_ratios);

    for (int i = 0; i < MAX_LEVELS; i++)
    {
        string keyname = get_alphan_key(isotope, target, i, database);
        string source_name = pwd + get_alphan_source_filename_from_key(keyname);
        if (!isRegularFile(source_name) || branch_ratios[i] == 0.)
        {
            continue;
        }
        get_neutron_energy_from_sundean(keyname, nEvents, neutron_energy_vals[i]);
    }

    TRandom3 rd(0);
    TRandom3 genor(0);
    gRandom->SetSeed(0);

    for (size_t i = 0; i < nEvents; i++)
    {
        int idx_decay = -1;
        double prob = rd.Uniform(0., 1.0);
        for (size_t j = 0; j < branch_ratios.size(); j++)
        {
            if (prob <= branch_ratios[j])
            {
                idx_decay = j;
                break;
            }
            else
            {
                prob -= branch_ratios[j];
            }
        }

        double cosTheta1 = genor.Uniform(-1., 1.);
        double azimuth1 = genor.Uniform(2.0 * M_PI);
        double theta1 = acos(cosTheta1);

        double Ek = neutron_energy_vals[idx_decay][i];
        double momentum1 = sqrt(pow(Ek / GeV + neutron_mass_c2 / GeV, 2) - pow(neutron_mass_c2 / GeV, 2));
        dirs p1(momentum1, theta1, azimuth1);

        vector<double> Eout;
        generate_deexcitation_gamma(key_mode[idx_decay], intensity[idx_decay], energy[idx_decay], Eout);
        int num_g = 1 + Eout.size();

        fprintf(stream, "%d\n", num_g);
        fprintf(stream, "1\t2112 0 0 %e %e %e %e 0\n",
                p1.x_norm, p1.y_norm, p1.z_norm, neutron_mass_c2 / GeV);

        for (int ig = 0; ig < num_g - 1; ig++)
        {
            cosTheta1 = genor.Uniform(-1., 1.);
            azimuth1 = genor.Uniform(2.0 * M_PI);
            theta1 = acos(cosTheta1);

            double Eg = Eout[ig];
            double momentum2 = Eg / GeV;
            dirs p2(momentum2, theta1, azimuth1);
            fprintf(stream, "1\t22 0 0 %e %e %e %e 0\n",
                    p2.x_norm, p2.y_norm, p2.z_norm, 0.);
        }
    }
}

// ==============================================
// Core Processing Functions
// ==============================================

/**
 * @brief Read and process NEUCBOT nuclear reaction code output
 */
int read_sundean_output(const char *input_filename, const char *output_filename)
{
    FILE *inputfile = fopen(input_filename, "r");
    vector<double> value_spec;

    if (strcmp(input_filename, "") != 0 && inputfile != nullptr)
    {
        char line[256];
        int is_print = 0;
        while (fgets(line, sizeof(line), inputfile))
        {
            string line_content = line;
            if (line_content[0] == '#')
                is_print++;
            if (is_print == 2 && line_content.size() > 1)
            {
                int spx = line_content.find_first_of(' ');
                double bin_content = atof(line_content.substr(spx + 1).data());
                value_spec.push_back(bin_content);
            }
        }
    }

    fclose(inputfile);
    string output_dir_path = output_filename;
    size_t found = output_dir_path.find_last_of("/\\"); // 查找最后一个 `/` 或 `\`
    string output_dir = "";
    if (found != std::string::npos)
    {
        output_dir = output_dir_path.substr(0, found); // 返回父目录路径
    }

    if (!directoryExists(output_dir))
    {
        createDirectory(output_dir);
    }

    TFile *outfile = new TFile(output_filename, "RECREATE");
    TH1D *hspec = new TH1D("hspec", "neutron spectrum", value_spec.size(), 0, value_spec.size() / 100);
    for (unsigned int i = 0; i < value_spec.size(); i++)
        hspec->SetBinContent(i + 1, value_spec[i]);
    hspec->Write();
    outfile->Close();
    delete outfile;

    return 0;
}

bool directoryExists(const std::string &path)
{
    struct stat info;
    if (stat(path.c_str(), &info) != 0)
    {
        return false; // 无法访问（可能不存在）
    }
    return (info.st_mode & S_IFDIR); // 如果是目录，返回 true
}

bool isRegularFile(const std::string &path)
{
    struct stat info;
    if (stat(path.c_str(), &info) != 0)
    {
        return false; // 无法访问（可能不存在）
    }
    return S_ISREG(info.st_mode); // 如果是普通文件，返回 true
}

bool createDirectory(const std::string &path)
{
    mode_t mode = 0755; // 权限：rwxr-xr-x
    if (mkdir(path.c_str(), mode) == 0)
    {
        return true; // 创建成功
    }
    if (errno == EEXIST)
    {
        return true; // 目录已存在
    }
    if (errno == ENOENT)
    { // 父目录不存在，尝试递归创建
        size_t pos = path.find_last_of("/\\");
        if (pos == std::string::npos)
        {
            return false;
        }
        std::string parent = path.substr(0, pos);
        if (!createDirectory(parent))
        {
            return false;
        }
        return (mkdir(path.c_str(), mode) == 0);
    }
    return false;
}

void TestDirectory(const string &source_name)
{
    size_t last_slash = source_name.find_last_of("/\\");
    if (last_slash != std::string::npos)
    {
        std::string dir_path = source_name.substr(0, last_slash);

        if (!directoryExists(dir_path))
        {
            if (!createDirectory(dir_path))
            {
                std::cerr << "Failed to create directory: " << dir_path << std::endl;
                // 处理错误（抛出异常、返回 false 等）
            }
            else
            {
                std::cout << "Created directory: " << dir_path << std::endl;
            }
        }
    }
}

/**
 * @brief Get neutron energies from processed NEUCBOT data
 */
void get_neutron_energy_from_sundean(const string &key,
                                     unsigned int nEvents, vector<double> &neutron_energy_vals)
{
    string source_name = pwd + get_alphan_source_filename_from_key(key);
    string data_name = pwd + get_alphan_root_filename_from_key(key);

    if (gSystem->AccessPathName(data_name.data()))
    {
        read_sundean_output(source_name.data(), data_name.data());
    }

    if (gSystem->AccessPathName(data_name.data()))
    {
        printf("Cannot create hist: %s\n", data_name.data());
        exit(1);
    }

    TFile *data_file = new TFile(data_name.data(), "READ");
    TH1D *hspec = (TH1D *)data_file->Get("hspec");
    for (unsigned int i = 0; i < nEvents; i++)
    {
        neutron_energy_vals.push_back(hspec->GetRandom());
    }

    delete data_file;
}

/**
 * @brief Get gamma transition data for single excited state
 */
void get_gamma_energy(const string &target, int state_idx,
                      vector<string> &key_mode, vector<double> &intensity, vector<double> &energy)
{
    string fname = get_alphan_target_filename(target, state_idx);
    FILE *inputfile = fopen(fname.c_str(), "r");

    if (fname != "" && inputfile != nullptr)
    {
        char line[256];
        int is_print = 0;
        while (fgets(line, sizeof(line), inputfile) != NULL)
        {
            string line_content = line;
            if (line_content.substr(0, 9) == "mother_id")
            {
                is_print = 1;
                continue;
            }
            if (is_print == 1 && line_content.size() > 1)
            {
                char *p = strtok(line, " \r\n");
                string content[3];
                int i = 0;
                while (p != NULL)
                {
                    content[i] = p;
                    i++;
                    p = strtok(NULL, " \n\r");
                }
                key_mode.push_back(content[0]);
                intensity.push_back(atof(content[1].c_str()));
                energy.push_back(atof(content[2].c_str()));
            }
        }
    }
}

/**
 * @brief Get gamma transition data for multiple excited states
 */
void get_gamma_energy(const string &target, int isotope, const string &database,
                      vector<string> key_mode[], vector<double> intensity[], vector<double> energy[],
                      vector<double> &branch_ratios)
{
    for (int state_idx = 0; state_idx < MAX_LEVELS; state_idx++)
    {
        string key = get_alphan_key(isotope, target, state_idx, database);
        string fname = pwd + get_alphan_source_filename_from_key(key);

        if (!isRegularFile(fname))
        {
            continue;
        }

        get_gamma_energy(target, state_idx,
                         key_mode[state_idx], intensity[state_idx], energy[state_idx]);

        FILE *inputfile = fopen(fname.c_str(), "r");
        string content[10];
        if (fname != "" && inputfile != nullptr)
        {
            char line[256];
            int is_print = 0;
            while (fgets(line, sizeof(line), inputfile) != NULL)
            {
                string line_content = line;
                if (line_content.size() > 1 && line_content.substr(0, 1) == "#" && is_print == 0)
                {
                    is_print = 1;
                    char *p = strtok(line, " \r\n");
                    int i = 0;
                    while (p != NULL)
                    {
                        content[i] = p;
                        i++;
                        p = strtok(NULL, " \r\n");
                    }
                }
            }
        }
        branch_ratios.push_back(atof(content[5].c_str()));
    }

    double sum_branch = accumulate(branch_ratios.begin(), branch_ratios.end(), 0.);
    for (size_t i = 0; i < branch_ratios.size(); i++)
    {
        branch_ratios[i] = branch_ratios[i] / sum_branch;
    }
}

/**
 * @brief Simulate gamma ray cascade de-excitation
 *
 * This is the core implementation of de-excitation via the mother id
 */
/**
 * @brief Simulate gamma ray cascade de-excitation process
 *
 * @param key_mode Vector of gamma transition identifiers (e.g. "0_1" for 0→1 transition)
 * @param intensity Vector of relative intensities for each transition
 * @param energy Vector of gamma energies (MeV) for each transition
 * @param Eout Output vector to store generated gamma energies
 *
 * @details
 * This function implements a Monte Carlo simulation of gamma cascades:
 * 1. Processes transitions in order of decreasing level numbers
 * 2. For each level, selects transitions based on:
 *    - Matching mother/daughter level IDs
 *    - Relative intensities (branching ratios)
 * 3. Continues until ground state is reached
 *
 * Algorithm:
 * 1. Scan transitions by level number (longest first)
 * 2. For each level, collect possible transitions
 * 3. Randomly select transition based on intensities
 * 4. Store gamma energy and continue to next level
 * 5. Repeat until ground state (level 0) is reached
 *
 * Note:
 * - key_mode format is "mother_daughter" (e.g. "2_0" for 2→0)
 * - Intensity values are normalized internally
 */
void generate_deexcitation_gamma(const vector<string> &key_mode,
                                 const vector<double> &intensity,
                                 const vector<double> &energy,
                                 vector<double> &Eout)
{
    const unsigned int len = energy.size();
    TRandom3 rd(0);
    vector<string> out_info;

    // length of mother id string
    for (unsigned int len_i = 2; len_i < 100; len_i++)
    {
        unsigned int num_level_smaller = 0;

        // this is a variable within one-level mother ids
        map<double, simpleObj> decay_candidates;

        // length of mother id's string vector
        for (unsigned int j = 0; j < len; j++)
        {
            // select enery level &&
            // (read first mother id || this mother id is following the last-level de-excitation)
            if (key_mode[j].size() == len_i &&
                (out_info.size() == 0 ||
                 (key_mode[j].find(out_info.back()) == 0 &&
                  key_mode[j].find(out_info.back()) != string::npos)))
            {
                simpleObj obj = {energy[j], key_mode[j]};
                // no de-excitation id recorded yet
                if (decay_candidates.size() == 0)
                {
                    decay_candidates[intensity[j]] = obj;
                }
                // already candidates recorded
                else
                {
                    map<double, simpleObj>::iterator iter = decay_candidates.end();
                    iter--;
                    // accumulate probability within one-level mother ids
                    double adding_prob = intensity[j] + (iter->first);
                    decay_candidates[adding_prob] = obj;
                }
            }
            if (key_mode[j].size() < len_i)
                num_level_smaller++;
        }

        if (num_level_smaller == len)
            break;
        if (decay_candidates.size() == 0)
            continue;

        size_t num_channels = decay_candidates.size();
        if (decay_candidates.rbegin()->first < 0.01)
        {
            double prob_uni = rd.Uniform(0., 1.0 * num_channels);
            int idx_decay = prob_uni;
            map<double, simpleObj>::iterator iter = decay_candidates.begin();
            for (int ii = 0; ii < idx_decay; ii++)
                iter++;

            out_info.push_back(iter->second.key_mode);
            Eout.push_back(iter->second.energy);
        }
        else
        {
            double prob = rd.Uniform(0., 1.);
            map<double, simpleObj>::iterator iter = decay_candidates.upper_bound(prob);

            out_info.push_back(iter->second.key_mode);
            Eout.push_back(iter->second.energy);
        }
    }
}

// ==============================================
// Utility Functions
// ==============================================

string get_alphan_source_filename_from_key(const string &key)
{
    string key_name = Form("%s.dat", key.c_str());
    return alpha_data_path + "sources/" + data_version_global + "/" + key_name;
}

string get_alphan_root_filename_from_key(const string &key)
{
    string key_name = Form("%s.root", key.c_str());
    string top_dir = alpha_data_path + "rootfile/" + data_version_global;
    TestDirectory(pwd + top_dir);
    return alpha_data_path + "rootfile/" + data_version_global + "/" + key_name;
}

string get_alphan_key(int isotope, const string &target,
                      int state_idx, const string &database)
{
    map<int, string> decayChainName;
    decayChainName[90232] = "Th232";
    decayChainName[92238] = "U238";
    decayChainName[92235] = "U235";
    decayChainName[82210] = "Pb210";
    decayChainName[89227] = "Ac227";
    decayChainName[1940] = "K40";

    if (target == "tot" || target == "all" || state_idx == -1)
    {
        // example: U238Alpha
        return Form("%sAlpha", decayChainName[isotope].c_str());
    }
    // example: b11_U238Alpha_state2_JENDL
    return Form("%s_%sAlpha_state%d_%s", target.c_str(), decayChainName[isotope].c_str(), state_idx, database.c_str());
}

string get_alphan_target_filename(const string &target, int state_idx)
{
    string key_name = Form("%s_state_%d.dat", target.c_str(), state_idx);
    return pwd + alpha_data_path + "sources/GammaInfoLib/" + key_name;
}
